## Daniel J. Mallinson and Peter K. Hatemi
## Reproduction Code for "The Effects of Information and Social Conformity on Opinion Change"
## 2017-08-17

rm(list = ls(all = TRUE))

install.packages("haven") #foreign package having trouble writing .dta file with strings
library(haven)

##############################################################################
###########################  Data Prep  ######################################
##############################################################################

experiment <- read.dta("raw_data.dta")

data <- experiment #Create working data object

#Fix incorrect variable name
library(plyr)
data <- rename(data, c(big5_neurotr="big5_neurot5r")) 

#Drop append, this is left over from STATA
#data$append <- NULL
myvar <- names(data) %in% c("append")
data <- data[!myvar]

#Strip out accessid
myvar <- names(data) %in% c("accessid")
data <- data[!myvar]

#Recode Revote so it is dichotomous and agreement is 1
data$revote <- mapvalues(data$revote, from= c(1,2), to= c(1,0))

#Make year of birth a numerical variable
data$yearofbirth <- as.numeric(data$yearofbirth)

#Create variable for age
data$age <- 2014-data$yearofbirth

#Recode variables that should be dichotomous
data$handed <- data$handed-1

##Create conservatism score

#Rename variables because of original naming error which mis-identified 
#flipped varaibles (resulting a liberalism score, not a conservatism score)
data <- rename(data, c(wpideo_1="wp1r"))
data <- rename(data, c(wpideo_2_r="wp2"))
data <- rename(data, c(wpideo_3_r="wp3"))
data <- rename(data, c(wpideo_4_r="wp4"))
data <- rename(data, c(wpideo_5_r="wp5"))
data <- rename(data, c(wpideo_6="wp6r"))
data <- rename(data, c(wpideo_7="wp7r"))
data <- rename(data, c(wpideo_8_r="wp8"))
data <- rename(data, c(wpideo_9="wp9r"))
data <- rename(data, c(wpideo_10="wp10r"))
data <- rename(data, c(wpideo_11_r="wp11"))
data <- rename(data, c(wpideo_12_r="wp12"))
data <- rename(data, c(wpideo_13_r="wp13"))
data <- rename(data, c(wpideo_14="wp14r"))
data <- rename(data, c(wpideo_15_r="wp15"))
data <- rename(data, c(wpideo_16="wp16r"))
data <- rename(data, c(wpideo_17="wp17r"))
data <- rename(data, c(wpideo_18="wp18r"))
data <- rename(data, c(wpideo_19_r="wp19"))
data <- rename(data, c(wpideo_20="wp20r"))
data <- rename(data, c(wpideo_21_r="wp21"))
data <- rename(data, c(wpideo_22="wp22r"))
data <- rename(data, c(wpideo_23_r="wp23"))
data <- rename(data, c(wpideo_24_r="wp24"))
data <- rename(data, c(wpideo_25_r="wp25"))
data <- rename(data, c(wpideo_26_r="wp26"))
data <- rename(data, c(wpideo_27="wp27r"))
data <- rename(data, c(wpideo_28_r="wp28"))
data <- rename(data, c(wpideo_29="wp29r"))
data <- rename(data, c(wpideo_30_r="wp30"))
data <- rename(data, c(wpideo_31_r="wp31"))
data <- rename(data, c(wpideo_32_r="wp32"))
data <- rename(data, c(wpideo_33="wp33r"))
data <- rename(data, c(wpideo_34_r="wp34"))
data <- rename(data, c(wpideo_35_r="wp35"))
data <- rename(data, c(wpideo_36_r="wp36"))
data <- rename(data, c(wpideo_37="wp37r"))
data <- rename(data, c(wpideo_38_r="wp38"))
data <- rename(data, c(wpideo_39="wp39r"))
data <- rename(data, c(wpideo_40_r="wp40"))
data <- rename(data, c(wpideo_41_r="wp41"))
data <- rename(data, c(wpideo_42_r="wp42"))
data <- rename(data, c(wpideo_43_r="wp43"))
data <- rename(data, c(wpideo_44="wp44r"))
data <- rename(data, c(wpideo_45="wp45r"))
data <- rename(data, c(wpideo_46="wp46r"))
data <- rename(data, c(wpideo_47="wp47r"))
data <- rename(data, c(wpideo_48="wp48r"))

#Flip the reverse coding
flip.cons <- c("wp1r", "wp6r", "wp7r", "wp9r", "wp10r", "wp14r", "wp16r", "wp17r", "wp18r", "wp20r", "wp22r", "wp27r", "wp29r", "wp33r", "wp37r", "wp39r", "wp44r", "wp45r", "wp46r", "wp47r", "wp48r")
data.flip.cons <- data[flip.cons]
for (i in 1:ncol(data.flip.cons)){
data.flip.cons[,i] <- mapvalues(data.flip.cons[,i], from= c(1,2), to=c(2,1))
} # flip the values
data[flip.cons] <- data.flip.cons[flip.cons] #replace in original data

data[,23:70] <- data[,23:70]-1 #Make a dichotomous indicator for each

data$conservatism <- rowSums(data[,23:70], na.rm=TRUE) #sum across the rows for each obs

## Create self-esteem score
flip.esteem <- c("esteem_3_r", "esteem_5_r", "esteem_8_r", "esteem_9_r", "esteem_10_r")
data.flip.esteem <- data[flip.esteem]
library(plyr)
for (i in 1:ncol(data.flip.esteem)){
  data.flip.esteem[,i] <- mapvalues(data.flip.esteem[,i], from= c(1,2,3,4,5), to=c(5,4,3,2,1))
} # flip the values
data[flip.esteem] <- data.flip.esteem[,1:5] #replace in original data

data[,116:125] <- data[,116:125]-1 #Make 0 to 4 scale
data$self_esteem <- rowSums(data[,116:125], na.rm=TRUE)

## Create Big 5 personality scores
flip.big5 <- c("big5_agree1r", "big5_extra2r", "big5_consc2r", "big5_neurot2r", "big5_agree3r", "big5_consc4r", "big5_extra5r", "big5_consc5r", "big5_neurot5r", "big5_agree6r", "big5_extra7r", "big5_neurot7r", "big5_open7r", "big5_agree8r", "big5_open9r")
data.flip.big5 <- data[flip.big5]
library(plyr)
for (i in 1:ncol(data.flip.big5)){
  data.flip.big5[,i] <- mapvalues(data.flip.big5[,i], from= c(1,2,3,4,5), to=c(5,4,3,2,1))
} # flip the values
data[flip.big5] <- data.flip.big5[,1:15] #Replace in original data

data[,72:115] <- data[,72:115]-1 #Make 0 to 4 scale

#Create additive scales for each of the Big 5 personality factors
data$extraversion <- rowSums(data[c("big5_extra1", "big5_extra2r", "big5_extra3", "big5_extra4", "big5_extra5r", "big5_extra6", "big5_extra7r", "big5_extra8")], na.rm=TRUE)
data$agreeableness <- rowSums(data[c("big5_agree1r", "big5_agree2", "big5_agree3r", "big5_agree4", "big5_agree5", "big5_agree6r", "big5_agree7", "big5_agree8r", "big5_agree9")], na.rm=TRUE)
data$conscientiousness <- rowSums(data[c("big5_consc1", "big5_consc2r", "big5_consc3", "big5_consc4r", "big5_consc5r", "big5_consc6", "big5_consc7", "big5_consc8", "big5_consc9")], na.rm=TRUE)
data$neuroticism <- rowSums(data[c("big5_neurot1", "big5_neurot2r", "big5_neurot3", "big5_neurot4", "big5_neurot5r", "big5_neurot6", "big5_neurot7r", "big5_neurot8")],na.rm=TRUE)
data$openness <- rowSums(data[c("big5_open1", "big5_open2", "big5_open3", "big5_open4", "big5_open5", "big5_open6", "big5_open7r", "big5_open8", "big5_open9r", "big5_open10")],na.rm=TRUE)

#Collapse initial survey Paterno vote into dichotomous variable
# 0 if disagree (4) or strongly disagree (5) with firing
# 1 if strongly agree (1), agree (2)
# 2 if neutral
data$paterno_recode <- 0
data$paterno_recode[data$paterno==1] <- 1
data$paterno_recode[data$paterno==2] <- 1
data$paterno_recode[data$paterno==3] <- 2

#Create recoded control group revote to reflect neutral category
data$revote_recode <- NA
data$revote_recode[data$revote==0 & data$revote_strength!=3] <- 0
data$revote_recode[data$revote==1 & data$revote_strength!=3] <- 1
data$revote_recode[data$revote_strength==3] <- 2
data[57,166] <- 2 #Consider missing revote_strength as neutral/don't know

#Create change variable
data$change <- 0
data$change[data$initialvote!=data$finalvote] <- 1
data$change[data$paterno_recode!=data$revote_recode] <- 1

#Create pro-change variable that groups neutrals with anti-firing respondents
data$pro_firing <- NA
data$pro_firing[data$paterno_recode==1] <- 1
data$pro_firing[data$paterno_recode==0] <- 0
data$pro_firing[data$paterno_recode==2] <- 0

##Add verbal_change back in to dataset
#This wasn't transfered to the original dataset
data$verbal_change <- NA

zero <- as.matrix(c(5,6, 10, 11, 13, 17, 20, 39, 40, 44, 46, 49, 57, 61, 62, 66, 67))
for(i in 1:nrow(zero)){
data$verbal_change[data$partnum==zero[i,1]] <- 0
}

one <- as.matrix(c(7, 8, 12, 16, 19, 21, 22, 41, 45, 48, 63, 64, 65))
for(i in 1:nrow(one)){
  data$verbal_change[data$partnum==one[i,1]] <- 1
}

### Strip out unnecessary X variables created by .csv saving
myvar <- names(data) %in% c("X", "X.1", "X.2")
data <- data[!myvar]


## Correct change category based on definition of change for control group ##
# Partnum 59, 25, 29, 36

data$change[data$partnum==59] <- 0
data$change[data$partnum==25] <- 0
data$change[data$partnum==29] <- 0
data$change[data$partnum==36] <- 0

write_dta(data, "reproduction_data.dta")

##############################################################################
###################  Analysis to Reproduce Manuscript ########################
##############################################################################

#install.packages("haven") #foreign package having trouble writing .dta file with strings
library(haven)

experiment <- read_dta("reproduction_data.dta")
data <- experiment

#Descriptive statistics for relevant variables
install.packages("psych")
library(psych)
describe(data[c("yearofbirth","age","male", "handed", "white", "maritalstatus", "schoolyear", "party", "ideology_self", "religiosity_self", "conservatism", "self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness")])
temp <- table(as.vector(data$maritalstatus))
names(temp)[temp == max(temp)]  # Mode for marital status
temp <- table(as.vector(data$schoolyear))
names(temp)[temp == max(temp)]  # Mode for school year
temp <- table(as.vector(data$party))
names(temp)[temp == max(temp)]  # Mode for political party
temp <- table(as.vector(data$religiosity_self))
names(temp)[temp == max(temp)]  # Mode for religiosity

##### Table of changes 
treat <- data[which(data$treatment==1),] # Pull out only treatment group
control <- data[which(data$treatment==0),] # Pull out only control group

#Pull only necessary variables
treat <- treat[c("partnum", "paterno", "initialvote", "verbal_agree_mode", "finalvote", "verbal_change", "change")]
control <- control[c("partnum", "paterno", "revote", "revote_strength")]

#write.csv(treat, file="treat_votes.csv")
#write.csv(control, file="control_votes.csv")

##### Comparison of time to complete control survey
control.time <- data[which(data$treatment==0),]
mean(control.time$complete_mins) # 24
mean(control.time[which(control.time$complete_mins < 60),]$complete_mins) #17
control.time.change <- control.time[which(control.time$change==1),]
mean(control.time.change$complete_mins)# 31
sd(control.time.change$complete_mins) # 31
control.time.nochange <- control.time[which(control.time$change==0),]
mean(control.time.nochange$complete_mins)# 23
sd(control.time.nochange$complete_mins) # 26
t.test(control.time.change$complete_mins, control.time.nochange$complete_mins) # p = 0.77

##Correlation between wilson-patterson ideology and self report ideology
library(Hmisc)
ideo.corr <- as.matrix(data[c("ideology_self", "conservatism")])
rcorr(ideo.corr)
cor(data$ideology_self, data$conservatism)

## Intercoder reliabilty for verbal change coding
library(reshape)
rater <- read.csv("verbal_coding.csv")

md <- melt(rater, id=c("partnum", "coder"))
new.rater <- cast(md, partnum~variable + coder)

library(irr)
irr <- kappam.fleiss(new.rater)
irr

#####################################################################################
############################## Table 1 ##############################################
#####################################################################################
#Logit model of treatment effect
treat.effect <- glm(change ~ treatment, data=data, family=binomial(link=logit))
summary(treat.effect)

#Odds ratio and CIs
ci1 <- confint(treat.effect, parm="treatment", method="Wald")
ctab1 <- cbind(est=coef(treat.effect)[2],ci1)
ortab1 <- exp(ctab1)
print(ortab1,digits=3)
100*(6.81-1) #581

# Robustness check with comparison of final vote and pre-test opinion for treatment group
data$change_recode <- data$change
data$change_recode[data$paterno_recode==2] <- NA
data$change_recode[data$paterno_recode!=2 & data$paterno_recode != data$finalvote] <- 1 #removes neutral (paterno_recode = 2)
data$change_recode[data$paterno_recode!=2 & data$paterno_recode == data$finalvote] <- 0 #removes neutral (paterno_recode = 2)

treat.robustness <- glm(change_recode ~ treatment, data=data, family=binomial(link=logit))
summary(treat.robustness)

ci2 <- confint(treat.robustness, parm="treatment", method="Wald")
ctab2 <- cbind(est=coef(treat.robustness)[2],ci2)
ortab2 <- exp(ctab2)
print(ortab2,digits=3)
100*(6.81-1) #581

#####################################################################################
######################## Figures S5a and S5b ########################################
#####################################################################################

## See "Votes for Both Conditions.xlsx" for the data supporting these figures

##### Bar chart of changes for treatment and control (SI Figures S2 and S3)
 
#1 = No Change (15/34)
#2 = Information Change Only (1/34)
#3 = Information Decision (5/34)
#4 = Information Decision and Discussion Change (5/34)
#5 = Discussion Change Only (8/34)

treat.table <- table(c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,2,3,3,3,3,3,4,4,4,4,4,5,5,5,5,5,5,5,5))
treat.table <- prop.table(treat.table)*100

jpeg("figures_s2.jpg", width=8, height=8, units="in", res=1000)
plot.new()
par(mai=c(1, 1,.5,.25))
colors <- c("#08519c", "#3182bd", "#6baed6", "#bdd7e7", "#eff3ff")
barplot(treat.table, ylab="Percentage", xlab="Treatment Condition", ylim=c(0,60),
	col=colors, names.arg="")
legend(4,55, legend=c("No Change", "Information Change", "Information Decision", "Information Decision and Discussion Change", "Discussion Change"), 
	cex=1, y.intersp=1, xjust=0.5, bty="n", fill=colors)
text(.75, 46, labels=c("44%"), pos=3)
text(1.95, 5, labels=c("3%"), pos=3)
text(3.1, 17, labels=c("15%"), pos=3)
text(4.3, 17, labels=c("15%"), pos=3)
text(5.5, 26, labels=c("24%"), pos=3)
dev.off()

#1 = No Change (13/24)
#2 = Information Change (1/24)
#3 = Information Decision (1/24)
#4 = Weaken Opinion (7/24)
#5 = Strengthen Opinion (2/24)
control.table <- table(c(1,1,1,1,1,1,1,1,1,1,1,1,1,2,3,4,4,4,4,4,4,4,5,5))
control.table <- prop.table(control.table)*100

jpeg("figure_s3.jpg", width=8, height=8, units="in", res=1000)
par(mai=c(1, 1,.5,.25))
colors <- c("#a50f15", "#de2d26", "#fb6a4a", "#fcae91", "#fee5d9")
barplot(control.table, ylab="Percentage", xlab="Control Condition", ylim=c(0,60),
	col=colors, names.arg="")
legend(4,55, legend=c("No Change", "Information Change", "Information Decision", "Weaken Opinion", "Strengthen Opinion"), 
	cex=1, y.intersp=1, xjust=0.5, bty="n", fill=colors)
text(.75, 56, labels=c("54%"), pos=3)
text(1.95, 6, labels=c("4%"), pos=3)
text(3.1, 6, labels=c("4%"), pos=3)
text(4.3, 31, labels=c("29%"), pos=3)
text(5.5, 10, labels=c("8%"), pos=3)
dev.off()

#####################################################################################
############################## Figure 3 #############################################
#####################################################################################

## See "Votes for Both Conditions.xlsx" for the data supporting these figures

#Visualize Overall Changes for Figure 3 based on collapsed categories from above
#1 Treatment Change 13/34 = 38 percent
#2 Treatment No Change 21/34 = 62 percent
#3 Control Change 2/24 = 8 percent
#4 Control No Change 22/24 = 92 percent

table.total <- matrix(c(22,21,2,13), ncol=2, byrow=TRUE)
colnames(table.total) <-  c("Control", "Treatment")
rownames(table.total) <-c("No Change", "Change")
table.total <- as.table(table.total)
table.total

chisq.test(table.total) # chi = 5.094, p = 0.024

tiff("fig3.tiff", width=8, height=8, units="in", res=500)
par(mai=c(1,1,.5,.25))
color <- c("blue3", "firebrick1") #color for paper and presentation
#color <- c("gray20", "gray80") #grayscale for pub
#color.legend <- c("gray80", "gray20")
x.label <- expression(paste(chi^2, "= 5.094, p < 0.05"))
barplot(round(prop.table(table.total, 2)*100, 0), ylab="Percentage", ylim=c(0,102), col=color, cex.axis=1.5, beside=TRUE, names.arg=c("Control", "Treatment"))
par(new=TRUE)
barplot(round(prop.table(table.total, 2)*100, 0), ylab="", ylim=c(0,102), col="black", cex.axis=1.5, beside=TRUE, names.arg=c("Control", "Treatment"), density=c(5), angle=c(45,135), axes=FALSE)
#legend(5,90, legend=c("Did Not Change Opinion", "Changed Opinion"), cex=1, y.intersp=1, xjust=0.5, bty="n", fill=color, density=c(5), angle=c(45,135))
title(xlab=x.label, cex.lab=1)
text(1.5, 92, labels=c("Change \n 92%"), pos=3)
text(2.5, 8, labels=c("No Change \n 8%"), pos=3)
text(4.5, 62, labels=c("Change \n 62%"), pos=3)
text(5.5, 38, labels=c("No Change \n 38%"), pos=3)
mtext(c("n=58"), side=1, adj=0, line=3)
dev.off()

#####################################################################################
############################## Table 2 ##############################################
#####################################################################################

#Create categories of behavior for anovas
treat$like_category <- NA
treat$like_category[treat$verbal_change > treat$change] <- 1
treat$like_category[treat$verbal_change < treat$change] <- 2
treat$like_category[treat$verbal_change==1 & treat$change==1] <- 3
treat$like_category[treat$verbal_change==0 & treat$change==0] <- 4

prop.table(table(treat$like_category))

#####################################################################################
############################## Table 3 ##############################################
#####################################################################################

## Pro versus anti-firing for entire set
#firing <- data[which(data$treatment==1),]
pro.firing <- data[which(data$paterno_recode==1),]
anti.firing <- data[which(data$paterno_recode==0),]
describe(pro.firing[c("age","male", "schoolyear", "conservatism", "ideology_self","self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness")])
describe(anti.firing[c("age","male", "schoolyear", "conservatism", "ideology_self","self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness")])
#t-tests of differences
t.test(pro.firing$age, anti.firing$age) #NS
t.test(pro.firing$schoolyear, anti.firing$schoolyear) #NS
t.test(pro.firing$conservatism, anti.firing$conservatism) #S.05
t.test(pro.firing$ideology_self, anti.firing$ideology_self) #NS
t.test(pro.firing$self_esteem, anti.firing$self_esteem) #NS
t.test(pro.firing$extraversion, anti.firing$extraversion) #NS
t.test(pro.firing$agreeableness, anti.firing$agreeableness) #NS
t.test(pro.firing$conscientiousness, anti.firing$conscientiousness) #NS
t.test(pro.firing$neuroticism, anti.firing$neuroticism)#NS
t.test(pro.firing$openness, anti.firing$openness) #NS

#chi-square for male
male.test <- rbind(pro.firing, anti.firing)
vote.table <- table(male.test$pro_firing, male.test$male)
prop.test(vote.table) #NS


#####################################################################################
############################## Table 4 ##############################################
#####################################################################################

## Differences between opinion changers and non-changers, overall
change <- data[which(data$change==1),]
no.change <- data[which(data$change==0),]

describe(change[c("age","male", "schoolyear", "confederates", "conservatism", "ideology_self","self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness", "pro_firing", "treatment")], skew=FALSE)
describe(no.change[c("age","male", "schoolyear", "confederates", "conservatism", "ideology_self","self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness", "pro_firing", "treatment")], skew=FALSE)
#t-tests of differences
t.test(change$age, no.change$age) #NS
t.test(change$schoolyear, no.change$schoolyear) #NS
t.test(change$confederates, no.change$confederates)#NS
t.test(change$conservatism, no.change$conservatism) #S.01
t.test(change$self_esteem, no.change$self_esteem) #NS
t.test(change$extraversion, no.change$extraversion) #NS
t.test(change$agreeableness, no.change$agreeableness) #NS
t.test(change$conscientiousness, no.change$conscientiousness) #.05
t.test(change$neuroticism, no.change$neuroticism)#S.05
t.test(change$openness, no.change$openness) #NS

#chi-square test for male
vote.table <- table(data$change, data$male)
prop.test(vote.table,correct=F) #NS


#####################################################################################
############################## Table 5 ##############################################
#####################################################################################

tchange <- data[which(data$change==1 & data$treatment==1),]
tno.change <- data[which(data$change==0 & data$treatment==1),]

describe(tchange[c("age","male", "schoolyear", "confederates", "conservatism", "self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness")], skew=FALSE)
describe(tno.change[c("age","male", "schoolyear", "confederates", "conservatism", "self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness")], skew=FALSE)
#t-tests of differences
t.test(tchange$age, tno.change$age) #NS
t.test(tchange$schoolyear, tno.change$schoolyear) #NS
t.test(tchange$confederates, tno.change$confederates)#NS
t.test(tchange$conservatism, tno.change$conservatism) #0.05
t.test(tchange$self_esteem, tno.change$self_esteem) #NS
t.test(tchange$extraversion, tno.change$extraversion) #NS
t.test(tchange$agreeableness, tno.change$agreeableness) #NS
t.test(tchange$conscientiousness, tno.change$conscientiousness) #.05
t.test(tchange$neuroticism, tno.change$neuroticism)#NS
t.test(tchange$openness, tno.change$openness) #NS

#chi-square test for male
vote.table <- table(data[which(data$treatment==1),]$change, data[which(data$treatment==1),]$male)
prop.test(vote.table,correct=F) #NS

#Robustness check - compare change and no change using collapsed pre-test opinion instead of verbal vote
t2change <- data[which(data$paterno_recode!=2 & data$treatment==1),]
t2no.change <- data[which(data$paterno_recode!=2 & data$treatment==1),]

t2change <- t2change[which(t2change$paterno_recode != t2change$finalvote),]
t2no.change <- t2no.change[which(t2no.change$paterno_recode == t2no.change$finalvote),]

describe(t2change[c("age","male", "schoolyear", "confederates", "conservatism", "self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness")], skew=FALSE)
describe(t2no.change[c("age","male", "schoolyear", "confederates", "conservatism", "self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness")], skew=FALSE)
#t-tests of differences
t.test(t2change$age, t2no.change$age) #NS
t.test(t2change$schoolyear, t2no.change$schoolyear) #NS
t.test(t2change$confederates, t2no.change$confederates)#0.05
t.test(t2change$conservatism, t2no.change$conservatism) #NS
t.test(t2change$self_esteem, t2no.change$self_esteem) #NS
t.test(t2change$extraversion, t2no.change$extraversion) #NS
t.test(t2change$agreeableness, t2no.change$agreeableness) #NS
t.test(t2change$conscientiousness, t2no.change$conscientiousness) #.10
t.test(t2change$neuroticism, t2no.change$neuroticism)#NS
t.test(t2change$openness, t2no.change$openness) #NS

#chi-square test for male
data.maletest <- data[which(data$paterno_recode!=2),]
data.maletest$pat_vote2 <- NA
data.maletest$pat_vote2[data.maletest$paterno_recode != data.maletest$finalvote] <- 1
data.maletest$pat_vote2[data.maletest$paterno_recode == data.maletest$finalvote] <- 0

vote.table <- table(data.maletest[which(data.maletest$treatment==1),]$pat_vote2, data.maletest[which(data.maletest$treatment==1),]$male)
prop.test(vote.table,correct=F) #NS

#####################################################################################
############################## Figure 4 #############################################
#####################################################################################

pro.firing <- data[which(data$pro_firing==1),]
anti.firing <- data[which(data$pro_firing==0),]

### Visualize distribution of conservatism for pro- and anti-paterno (Figure 4)
density.profiring <- density(pro.firing$conservatism)
density.antifiring <- density(anti.firing$conservatism)
density.overall <- density(data$conservatism)

tiff("fig4.tiff", height=6, width=6, units="in", res=500)
par(mai=c(1, 1,.3,.25))
plot(density.profiring, col="blue3", xlab="Kernel Densities of Conservatism for Pro- and Anti-Firing Groups", main="", ylim=c(0,0.07), xlim=c(-10,50), axes=T)
polygon(density.profiring, col=adjustcolor("black",alpha.f=0.7), density=20, angle=45)
polygon(density.profiring, col=adjustcolor("blue3",alpha.f=0.7))
lines(density.antifiring, col="firebrick1", ylim=c(0,3), xlim=c(-10, 50))
polygon(density.antifiring, col=adjustcolor("black",alpha.f=0.7), density=20, angle=75)
polygon(density.antifiring, col=adjustcolor("firebrick1",alpha.f=0.7))
text(18, 0.05, "Pro-Firing", col="blue4")
text(30, 0.03, "Anti-Firing", col="firebrick1")
dev.off()

#####################################################################################
##############################    Table S3    #######################################
#####################################################################################

treat.not_random <- treatment.all[1:16,]
treat.random <- treatment.all[17:34,]
#Table of descriptives
describe(treat.not_random[c("age","male", "schoolyear", "conservatism", "self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness", "pro_firing")])
describe(treat.random[c("age","male", "schoolyear", "conservatism", "self_esteem", "extraversion", "agreeableness", "conscientiousness", "neuroticism", "openness", "pro_firing")])

#Testing differences between groups
t.test(treat.not_random$age, treat.random$age) #NS
t.test(treat.not_random$schoolyear, treat.random$schoolyear) #NS
t.test(treat.not_random$conservatism, treat.random$conservatism) #NS
t.test(treat.not_random$self_esteem, treat.random$self_esteem) #NS
t.test(treat.not_random$extraversion, treat.random$extraversion) #NS
t.test(treat.not_random$agreeableness, treat.random$agreeableness) #NS
t.test(treat.not_random$conscientiousness, treat.random$conscientiousness) #NS
t.test(treat.not_random$neuroticism, treat.random$neuroticism)#S.05
t.test(treat.not_random$openness, treat.random$openness) #NS

#chi-square tests for male, and paterno
treatment.all$nonrandom <- 1
treatment.all[17:34, 171] <- 0
vote.table <- table(treatment.all$male, treatment.all$nonrandom)
prop.test(vote.table,correct=F) #S.05
vote.table <- table(treatment.all$pro_firing, treatment.all$nonrandom)
prop.test(vote.table, correct=F) #NS

#####################################################################################
##############################    Figure S3    ######################################
#####################################################################################

data.random <- data[17:58,]

#Visualize Overall Changes for non-random data
#1 Treatment Change 6/18 = 33 percent
#2 Treatment No Change 12/18 = 67 percent
#3 Control Change 2/22 = 8 percent
#4 Control No Change 22/24 = 92 percent

rand.total <- matrix(c(22,12,2,6), ncol=2, byrow=TRUE)
colnames(rand.total) <-  c("Control", "Treatment")
rownames(rand.total) <-c("No Change", "Change")
rand.total <- as.table(rand.total)
rand.total

chisq.test(rand.total) #x = 2.71, p = 0.1

jpeg("change_barplot_total_random.jpg", width=8, height=8, units="in", res=1000)
par(mai=c(1,1,.5,.25))
color <- c("royalblue4", "firebrick") #color for paper and presentation
#color <- c("gray20", "gray80") #grayscale for pub
#color.legend <- c("gray80", "gray20")
x.label <- expression(paste(chi^2, "= 2.705, p = 0.10"))
barplot(round(prop.table(rand.total, 2)*100, 0), ylab="Percentage", ylim=c(0,100), col=color, cex.axis=1.5, beside=TRUE, names.arg=c("Control", "Treatment"))
legend(5,95, legend=c("Did Not Change Opinion", "Changed Opinion"), cex=1, y.intersp=1, xjust=0.5, bty="n", fill=color)
title(xlab=x.label, cex.lab=1)
text(1.5, 94, labels=c("92%"), pos=3)
text(2.5, 10, labels=c("8%"), pos=3)
text(4.5, 69, labels=c("67%"), pos=3)
text(5.5, 35, labels=c("33%"), pos=3)
mtext(c("n=42"), side=1, adj=0, line=3)
dev.off()
